/*
 * MAEquations.h by Michel J. Anders (c)2013
 * from https://github.com/sambler/osl-shaders
 *
 * license: cc-by-sa
 *
 * original script from -
 * http://blenderthings.blogspot.nl/p/a-small-osl-libray.html
 *
 */

// cubic roots adapted for OSL from http://van-der-waals.pc.uni-koeln.de/quartic/quartic.html
// so now we have a translation from Fortran -> C -> OSL. It doesn't look that well, but it works :-)

float CBRT(float Z) { return abs(pow(abs(Z),1.0/3.0)) * sign(Z); }

/*-------------------- Global Function Description Block ----------------------
*
*     ***CUBIC************************************************08.11.1986
*     Solution of a cubic equation
*     Equations of lesser degree are solved by the appropriate formulas.
*     The solutions are arranged in ascending order.
*     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
*     ******************************************************************
*     A(0:3)      (i)  vector containing the polynomial coefficients
*     X(1:L)      (o)  results
*     L           (o)  number of valid solutions (beginning with X(1))
*     ==================================================================
*   17-Oct-2004 / Raoul Rausch
*       Conversion from Fortran to C
*       09-Jan-2012 / Michel Anders (varkenvarken)
*               Conversion from C to OSL
*
*-----------------------------------------------------------------------------
*/

int cubic(float A[4], float X[3], int L)
{
    float PI = 3.1415926535897932;
    float THIRD = 1./3.;
    float U[3],W, P, Q, DIS, PHI;
    int i;

    // ====determine the degree of the polynomial ====

    if (A[3] != 0.0)
    {
        //cubic problem
        W = A[2]/A[3]*THIRD;
        P = pow((A[1]/A[3]*THIRD - pow(W,2)),3);
        Q = -.5*(2.0*pow(W,3)-(A[1]*W-A[0])/A[3] );
        DIS = pow(Q,2)+P;
        if ( DIS < 0.0 )
        {
            //three real solutions!
            //Confine the argument of ACOS to the interval [-1;1]!
            PHI = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
            P=2.0*pow((-P),(5.e-1*THIRD));
            for (i=0;i<3;i++)
                U[i] = P*cos((PHI+2*((float)i)*PI)*THIRD)-W;
            X[0] = min(U[0], min(U[1], U[2]));
            X[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
            X[2] = max(U[0], max(U[1], U[2]));
            L = 3;
        }
        else
        {
            // only one real solution!
            DIS = sqrt(DIS);
            X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
            L=1;
        }
    }
    else if (A[2] != 0.0)
    {
        // quadratic problem
        P = 0.5*A[1]/A[2];
        DIS = pow(P,2)-A[0]/A[2];
        if (DIS > 0.0)
        {
            // 2 real solutions
            X[0] = -P - sqrt(DIS);
            X[1] = -P + sqrt(DIS);
            L=2;
        }
        else
        {
            // no real solution
            L=0;
        }
    }
    else if (A[1] != 0.0)
    {
        //linear equation
        X[0] =A[0]/A[1];
        L=1;
    }
    else
    {
        //no equation
        L=0;
    }
/*
*     ==== perform one step of a newton iteration in order to minimize
*          round-off errors ====
*/
    for (i=0;i < L;i++)
    {
        X[i] = X[i] - (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3])))
            /(A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3]));
    }

    return 0;
}


point cubicspline(float t, point P0, point P1, point P2, point P3)
{
    return (1-t)*(1-t)*(1-t)*P0
            + 3*(1-t)*(1-t)*t*P1
            + 3*(1-t)*t*t*P2
            + t*t*t*P3;

    // how to group all factors of the powers of t
    //
    // (1-2t+t2-t+2t2+t3)P0 => (1-3t+3t2+t3)P0
    // (3t-6t2+3t3)P1
    // (3t2-3t3)P2
    // t3P3
    //
    // 1 * ( P0 )
    // t * ( 3P0+3P1 )
    // t2* ( 3P0-6P1+3P2 )
    // t3* ( P0+3P1-3P2+P3)
}

// calculate the closest distance from Pos to a quadratic bezier curve
// the curve is defined by the 3 points P0, P1 and P2

int splinedist(point p0, point p1, point p2, point Pos, float d, float tc)
{

    point P0 = p0;
    point P1 = p1;
    point P2 = p2;

    // following definitions are for the four polynomic coefficients for the well known
    // equation dB/dt . (Pos-B)  (i.e. the inproduct of the tangent to the bezier and the
    // difference vector from the point under considertion to the Bezier curve.
    // If the difference vector is perpendicular to the tangent we have found a closest point
    // on the Bezier curve.
    // The stuff below is generated by a script and no effort is spent on collecting factors.
    // We let the OSL compiler worry about that :-)

    float t0 = -2*P0[0]*P0[0]+-2*P1[0]*Pos[0]+2*P0[0]*P1[0]+2*P0[0]*Pos[0]
                +-2*P0[1]*P0[1]+-2*P1[1]*Pos[1]+2*P0[1]*P1[1]+2*P0[1]*Pos[1]
                +-2*P0[2]*P0[2]+-2*P1[2]*Pos[2]+2*P0[2]*P1[2]+2*P0[2]*Pos[2];

    float t1 = -4*P0[0]*P1[0]+-4*P0[0]*P1[0]+-4*P0[0]*P1[0]+-2*P0[0]*Pos[0]
                +-2*P2[0]*Pos[0]+2*P0[0]*P0[0]+2*P0[0]*P2[0]+4*P0[0]*P0[0]
                +4*P1[0]*P1[0]+4*P1[0]*Pos[0]+-4*P0[1]*P1[1]+-4*P0[1]*P1[1]
                +-4*P0[1]*P1[1]+-2*P0[1]*Pos[1]+-2*P2[1]*Pos[1]+2*P0[1]*P0[1]
                +2*P0[1]*P2[1]+4*P0[1]*P0[1]+4*P1[1]*P1[1]+4*P1[1]*Pos[1]
                +-4*P0[2]*P1[2]+-4*P0[2]*P1[2]+-4*P0[2]*P1[2]+-2*P0[2]*Pos[2]
                +-2*P2[2]*Pos[2]+2*P0[2]*P0[2]+2*P0[2]*P2[2]+4*P0[2]*P0[2]
                +4*P1[2]*P1[2]+4*P1[2]*Pos[2];

    float t2 = -8*P1[0]*P1[0]+-4*P0[0]*P0[0]+-4*P0[0]*P2[0]+-4*P1[0]*P1[0]
                +-2*P0[0]*P0[0]+-2*P0[0]*P2[0]+2*P0[0]*P1[0]+2*P1[0]*P2[0]
                +4*P0[0]*P1[0]+4*P0[0]*P1[0]+4*P1[0]*P2[0]+8*P0[0]*P1[0]
                +-8*P1[1]*P1[1]+-4*P0[1]*P0[1]+-4*P0[1]*P2[1]+-4*P1[1]*P1[1]
                +-2*P0[1]*P0[1]+-2*P0[1]*P2[1]+2*P0[1]*P1[1]+2*P1[1]*P2[1]
                +4*P0[1]*P1[1]+4*P0[1]*P1[1]+4*P1[1]*P2[1]+8*P0[1]*P1[1]
                +-8*P1[2]*P1[2]+-4*P0[2]*P0[2]+-4*P0[2]*P2[2]+-4*P1[2]*P1[2]
                +-2*P0[2]*P0[2]+-2*P0[2]*P2[2]+2*P0[2]*P1[2]+2*P1[2]*P2[2]
                +4*P0[2]*P1[2]+4*P0[2]*P1[2]+4*P1[2]*P2[2]+8*P0[2]*P1[2];

    float t3 = -4*P0[0]*P1[0]+-4*P0[0]*P1[0]+-4*P1[0]*P2[0]+-4*P1[0]*P2[0]
                +2*P0[0]*P0[0]+2*P0[0]*P2[0]+2*P0[0]*P2[0]+2*P2[0]*P2[0]
                +8*P1[0]*P1[0]+-4*P0[1]*P1[1]+-4*P0[1]*P1[1]+-4*P1[1]*P2[1]
                +-4*P1[1]*P2[1]+2*P0[1]*P0[1]+2*P0[1]*P2[1]+2*P0[1]*P2[1]
                +2*P2[1]*P2[1]+8*P1[1]*P1[1]+-4*P0[2]*P1[2]+-4*P0[2]*P1[2]
                +-4*P1[2]*P2[2]+-4*P1[2]*P2[2]+2*P0[2]*P0[2]+2*P0[2]*P2[2]
                +2*P0[2]*P2[2]+2*P2[2]*P2[2]+8*P1[2]*P1[2];

    float A[4] = {t0,t1,t2,t3};
    float T[3] ;
    int n ;
    cubic(A,T,n);
    d = 1e6;
    // cubic() will return 0 , 1 or 3 values for t
    // we are only interested in values that lie in the interval [0,1]
    // for those we calculate the position on the curve and check whether
    // we have found the shortest distance.
    int found = 0;
    while(n>0){
        n--;
        if(T[n]>=0 && T[n]<=1){
            float t = T[n];
            found = 1;
            float dd = distance((1-t)*(1-t)*P0 + 2*(1-t)*t *P1 + t*t*P2, Pos);
            if (dd < d) {
                d = dd;
                tc = t;
            }
        }
    }
    return found;
}


